arXiv:1508.05937vl [cond-mat.other] 24 Aug 2015 


A Topological Framework for Local Structure Analysis in Condensed Matter 


Emanuel A. Lazar^, Jian Han^, David J. Srolovitz^’^ 

^Department of Materials Seienee and Engineering 
^Department of Meehanieal Engineering and Applied Meehanies 
University of Pennsylvania, Philadelphia, Pennsylvania 191 Of, USA. 

(Dated: August 26, 2015) 

Physical systems are frequently modeled as sets of points in space, each representing the position 
of an atom, molecule, or mesoscale particle. As many properties of such systems depend on the 
underlying ordering of their constituent particles, understanding that structure is a primary objective 
of condensed matter research. Although perfect crystals are fully described by a set of translation 
and basis vectors, real-world materials are never perfect, as thermal vibrations and defects introduce 
significant deviation from ideal order. Meanwhile, liquids and glasses present yet more complexity. 
A complete understanding of structure thus remains a central, open problem. Here we propose 
a unified mathematical framework, based on the topology of the Voronoi cell of a particle, for 
classifying local structure in ordered and disordered systems that is powerful and practical. We 
explain the underlying reason why this topological description of local structure is better suited for 
structural analysis than continuous descriptions. We demonstrate the connection of this approach 
to the behavior of physical systems and explore how crystalline structure is compromised at elevated 
temperatures. We also illustrate potential applications to identifying defects in plastically deformed 
polycrystals at high temperatures, automating analysis of complex structures, and characterizing 
general disordered systems. 


Condensed matter systems are often abstracted as 
large sets of points in space, each representing the po¬ 
sition of an atom, molecule, or mesoscale particle. Two 
challenges frequently encountered when studying systems 
at this scale are classifying and identifying local struc¬ 
ture. Simulation studies of nucleation, crystallization, 
and melting, for example, as well as those of defect mi¬ 
gration and transformation, require a precise understand¬ 
ing of which particles are associated with which phases, 
and which are associated with defects. As these systems 
are abstracted as large point sets, these dual challenges 
of classifying and identifying structure reduce to ones of 
understanding arrangements of points in space. 

A primary difficulty in classifying structure in spatial 
point sets arises from a tension between a desire for com¬ 
pleteness and the necessity for practicality. The local 
neighborhood of a particle within an ensemble of par¬ 
ticles can be completely described by a list of relative 
positions of each of its neighbors. However, while such 
a list of coordinates is complete in some sense, this raw 
data provides little direct insight, leaving us wanting for a 
practical and more illuminating description. This tension 
is often mediated by the choice of an “order parameter”, 
which distills structural data into a single number or vec¬ 
tor, and which is constructed to be both informative and 
computationally tractable [1]. 

A central limitation of conventional order parameters 
is exhibited in degeneracies that arise in describing neigh¬ 
borhoods that are structurally distinct but which map 
to identical order-parameter values. Some order parame¬ 
ters classify particles in face-centered cubic (FCC) and 
body-centered cubic (BCC) crystals identically, while 
others classify particles in FCC and hexagonal close- 
packed (HCP) crystals identically [2]. Similarly, parti¬ 
cles located near defects in a low-temperature crystal can 


have order-parameter values identical to those of parti¬ 
cles in a high-temperature defect-free crystal. These de¬ 
generacies point to an inherent incompleteness in such 
order-parameter classifications of local structure. Conse¬ 
quently, different order parameters are necessary to study 
different systems BE]. 



FIG. 1. The frame of a Voronoi cell of a central particle (blue), 
surrounded by its nearest neighbors (gold). The topology of 
the Voronoi cell captures structural information about the 
local neighborhood. 

In this paper we propose a mathematical framework 
to classify local structure that avoids much of the degen¬ 
eracy encountered in other approaches and which, there¬ 
fore, is equally applicable to all ordered and disordered 
systems of particles. More specifically, the local struc¬ 
ture around a particle is classified using the topology of 
its Voronoi cell (see Fig. [^. Families of Voronoi cell 
topologies are constructed by considering those topolo¬ 
gies that can result from infinitesimal perturbations of 
an ideal structure. We demonstrate that this classifica¬ 
tion scheme is consistent with the manner in which local 
ordering changes in high-temperature single crystals as 
the temperature is raised toward their melting points. 
We highlight a potential use of this approach for visu¬ 
alizing defects in crystalline solids at high temperatures. 
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and contrast it with previous methods. We then demon¬ 
strate an application of this approach to the automated 
analysis of the evolution of complex structures, where 
conventional methods are often inadequate. Finally, we 
show an application in which this approach is used to 
provide robust statistical-structural descriptors for char¬ 
acterizing disordered systems. 


I. THEORY 

A. The Configuration Space of Local Structure 

A deeper understanding of local structure can be devel¬ 
oped through consideration of all possible arrangements 
of neighbors of a central particle. The local neighborhood 
of a particle within an ensemble of particles can be com¬ 
pletely described by a vector of relative positions of its 
n nearest neighbors: x = (ri, r2,..., r^), where is the 
relative position of the ith neighbor of a central particle. 
For suitably large n, any question about the local neigh¬ 
borhood of a particle can be answered through complete 
knowledge of x. We use C{n) to denote the configuration 
space of all possible arrangements of n nearest neighbors: 

C(n) = {(ri,r 2 ,...,r„) : Fj e (1) 

Each point in C{n) thus corresponds to a specific local ar¬ 
rangements of particles. Figure |^d) provides a schematic 
of C(n) and highlights points corresponding to local ar¬ 
rangements of BCC, FCC, HCP, and diamond structures. 
As defined in Eq.[^ the dimension of C(n) is 3n; ignoring 
rotations and scaling reduces the dimension of C(n) by 4. 
Ignoring permutations of the n neighbors and disallowing 
multiplicities further changes the geometry and topology 
of C(n), but not its dimension. 

Order parameters can be thought of as functions that 
map C{n) to a lower-dimensional order-parameter space; 
order-parameter spaces most commonly used are 
where d is substantially smaller than 3n. Each choice 
of order parameter results in a different subdivision of 
C{n) into regions on which that order parameter is con¬ 
stant; for real-valued continuous functions, these regions 
are commonly known as level sets. To help understand 
the degeneracy observed in continuous order-parameter 
methods, consider that for every continuous mapping 0 
from an unbounded high-dimensional space to a lower¬ 
dimensional space, there exist points xi,X 2 arbitrarily 
far apart, but for which (/)(xi) and 0 (x 2 ) are identical 
[3]. The continuity of an order parameter thus entails 
the kind of degeneracy highlighted above. In contrast, 
discrete order parameters are not subject to this limi¬ 
tation, as distances between points with identical order- 
parameter values can be bounded. This motivates the 
question of how to reasonably subdivide C{n). We now 
show that Voronoi topology offers one such approach. 




FIG. 2. (a-c) Voronoi cells of particles in BCC, FCC, and 
HCP crystals. Vertices at which more than four Voronoi cells 
meet are marked by red circles. Small perturbations of the 
particle positions result in topological changes near these ver¬ 
tices. (d) Schematic of C(n), the space of all possible configu¬ 
rations of n neighbors. This space can be divided into regions 
on which the Voronoi cell topology is constant. The topology 
of a point that lies on the boundaries of multiple regions is 
unstable, and infinitesimal perturbations of the neighbors will 
result in a change of topology. The inset shows the neighbor¬ 
hood around xpcc; unshaded regions indicate primary types, 
while shaded ones indicate secondary types. 


B. Voronoi Topology 

For a fixed set of particles, the Voronoi cell of a central 
particle is the region of real space closer to that parti¬ 
cle than to any other [4]. Figure illustrates a central 
particle, its Voronoi cell, and fifteen neighboring parti¬ 
cles. Two particles whose Voronoi cells share a face are 
called neighbors. We identify two Voronoi cells as having 
the same topology, or topological type, if there exists 
a one-to-one correspondence between their sets of faces 
that preserves adjacency. 

The topology of the Voronoi cell of a particle describes 
the manner in which neighbors of a particle are arranged 
relative to it. An n-sided face, for example, indicates a 
pair of particles which have n neighbors in common. The 
topology of a Voronoi cell thus provides a robust descrip¬ 
tion of how neighbors are arranged not only relative to a 
central particle, but also to one another. In this sense, it 
is a good description of local structure. 

Voronoi cell topology also provides a natural decom¬ 
position of C{n) into regions in which the Voronoi cell 
topology is constant, as illustrated in Fig. |^d). We con¬ 
sider this decomposition natural because it allows us to 
coarse-grain the effects of small perturbations on local 
structure. Small perturbations of the particle coordi¬ 
nates correspond to small displacements in C(n), and 
since the Voronoi topology does not change for almost all 
points under small perturbations, these small perturba¬ 
tions, which are often unimportant, are naturally ignored 
without the introduction of an artificial cut-off [2]. 
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Voronoi cell topology was first introduced by Bernal 
and others to study the atomic structure of liquids ISHZl, 
and has been subsequently applied to study a wide range 
of condensed matter systems, including random sphere 
packings miH], finite-temperature crystals [9] , and metal¬ 
lic glasses m- In those studies, however, the topology 
of a cell was characterized by counting its types of faces 
(e.g., triangles and quadrilaterals), though it ignored the 
way in which those faces are arranged. While this limited 
description has been used to study some aspects of crys¬ 
tallization m , it cannot distinguish particles whose local 
environments are FCC from those whose local environ¬ 
ments are HCP, as both Voronoi cells have twelve four¬ 
sided faces. In previous work [Tans], the authors have 
shown how to use a graph-tracing algorithm introduced 
by Weinberg la to efficiently compute strings which en¬ 
code a complete description of the Voronoi cell topology; 
see Methods for further details. 

A second limitation arising in traditional Voronoi ap¬ 
proaches results from abrupt changes in topology due to 
small geometric perturbations. Consider, for example, 
that Voronoi cells of particles in FCC and HCP crys¬ 
tals are topologically unstable - since some vertices are 
shared by more than four Voronoi cells (see Fig. §, in- 
finitesimal perturbations of the particle positions, such 
as those arising from non-hydrostatic strain or thermal 
vibrations, will change their topology [15]. This prob¬ 
lem has been sufficiently challenging to limit the util¬ 
ity of conventional Voronoi approaches in studying even 
slightly perturbed crystal structures [2]. This problem 
can be solved through the classification of topological 
types described in the following section. 


C. Theory of A-types 

In this section we show how topological types can be 
classified using the approach developed in the previous 
two sections. On a basic level, every arrangement of 
neighbors relative to a central particle can be described 
by its Voronoi cell topology. Families of topological types 
associated with a particular structure can then be defined 
as sets of types obtained through infinitesimal perturba¬ 
tions of that structure. This classification scheme enables 
a description of the effects of small strains and thermal vi¬ 
brations on local structure, and provides a robust frame¬ 
work suitable for theoretical and numerical analysis. 

Every local arrangement of neighbors A is described by 
a distinct point xa in C(n), and subsequently corresponds 
to a unique Voronoi cell topology E[xa]. For example, 
if A = BCC, then xa = xbcc describes a particle that 
has the same local environment as a particle in a perfect 
BCC crystal; its Voronoi cell topology E[xbcc] is the 
truncated octahedron, illustrated in Fig. |^a). 

A suitable distance function on C{n) allows us to de¬ 
fine sets of topological types associated with infinitesimal 
perturbations as follows. We let be a ball of ra¬ 

dius e centered at x. This region of C{n) corresponds to 


configurations obtained through small perturbations of a 
particle and its neighbors, where e controls the magni¬ 
tude of such a perturbation. The set of topological types 
obtained from all possible perturbations of x with mag¬ 
nitude smaller than e is denoted V[B^{x.)]. We define the 
family of topological types associated with infinitesimal 
perturbations of A as the limiting set: 


( 2 ) 


In more physical terms, Bx is the set of all topological 
types that can be obtained through arbitrarily small per¬ 
turbations of a central particle and its neighbors. The 
Voronoi cell topology of points in the interior of a re¬ 
gion in C{n) remains unchanged by small perturbations. 
In contrast, points such as xpcc and xhcp are located 
at the boundaries of multiple regions, and small pertur¬ 
bations entail a change in Voronoi cell topology. Thus, 
J^FCC and Jmcp consist of multiple topological types, 
whereas Bbcc^ located in the interior of a region, con¬ 
sists of a single type. If a topological type is in Bx, then 
we say that it is a A-type. Note that a topological type 
can belong to multiple families; this indeterminacy will 
be considered below. This classification of A-types allows 
us to account for topological instability without modify¬ 
ing the sample data by collapsing edges or faces using ad 
hoc criteria (e.g., cut-offs) [2j[9]. 

Among A-types, a further distinction can be drawn 
based on the manner in which the Voronoi cell topol¬ 
ogy subdivides C{n). Using a suitable volume measure 
Vol, we define the ideal frequency /a(p) of a topological 
type r relative to xa as follows: 


fx{T) := lim 

e^O 


Vol (y-i[r]nBAxA)) 
Vol(B,(x;,)) 


(3) 


where V~^[t] is the set of points in C{n) whose Voronoi 
cell have topology r. If /a(p) > 0, we call r a primary 
A-type; if fx{^) = 0, we call it a secondary A-type. The 
inset in Fig. ^d) highlights a number of regions incident 
with xpcc- Some of those regions meet xpcc at finite 
solid angles; therefore, their fractional volumes within 
an e-ball converge to positive values as e ^ 0; these are 
primary FCC-types. In contrast, fractional volumes tend 
to zero as e ^ 0 for other regions which meet xpcc at 
cusps; these are secondary FCC-types. 

The distinction between primary and secondary types 
appears to result from the manner in which topological 
instabilities resolve when perturbed. To illustrate this 
distinction. Fig. [^a) shows an unstable vertex shared by 
six Voronoi cells in FCC or HCP crystals; such vertices 
are marked by red circles in Figs. I^b) and (c). Fig¬ 
ure |^b) depicts the most common manner in which such 
an unstable vertex resolves upon small perturbations of 
neighboring particles m- In this resolution, a new four¬ 
sided face is formed between two non-adjacent Voronoi 
cells; all unstable vertices resolve in this manner in pri¬ 
mary types. A less common resolution can also occur as 
a result of cooperative motion of neighboring particles. 



4 



FIG. 3. (a) An unstable vertex shared by six Voronoi cells in 
FCC or HOP crystals; such vertices are marked by red circles 
in Figs. I^b) and (c). A small perturbation will cause the 
vertex to resolve into either (b) a four-sided face, or (c) a pair 
of adjacent triangular faces; these resolutions are associated 
with primary and secondary types, respectively. 

In this resolution, depicted in Fig. [^c), two triangular 
faces are created m; secondary types can include such 
resolutions. 

Determining J^x is feasible through consideration of 
all possible ways in which unstable vertices can resolve. 
For example, the ideal FCC Voronoi cell, illustrated in 
Fig. ib), has six unstable vertices. In primary types, 
each such vertex resolves in a manner illustrated in 
Fig. ib), in one of three directions. More specifically, 
the unstable vertex can transform in such a way that the 
central cell gains a square face, or else gains an edge in 
one of two directions. We consider all possible combina¬ 
tions of these resolutions over the six unstable vertices, 
and calculate the topological types of the resulting poly- 
hedra using the algorithm developed in [13]; a total of 
44 distinct topological types occur in this manner. In 
secondary types, unstable vertices can also resolve in the 
manner illustrated in Fig. |^c), or else remain unstable. 
An additional 6250 topological types can occur in this 
manner. A similar approach can be followed to deter¬ 
mine Tx for other structures. Additional details can be 
found in the supplementary material. 


II. FINITE-TEMPERATURE CRYSTALS 

The proposed distinction between primary and sec¬ 
ondary types is supported by atomistic simulation. We 
studied the atomic structure of three model materials, 
BCC tungsten [16], FCC copper [17], and HCP mag¬ 
nesium m, at elevated temperatures using molecular 
dynamics (MD) in the NPT ensemble [19]. Simulated 
systems contained 1,024,000, 1,372,000, and 1,029,600 
atoms, respectively, in periodic supercells. In each sim¬ 
ulation, a defect-free crystal was heated from T = 0 in 
increments of 50 K and equilibrated for 50 ps at each 
temperature. Figure shows how the distribution of 
topological types changes with temperature; each curve 
indicates the frequency of a single topological type. 

Types can be grouped according to the shape of their 
frequency curves. Frequencies of one group of types ap¬ 
proach finite values as T ^0, and change very little with 
temperature (blue curves in Fig.|^. Frequencies of a sec¬ 
ond group rapidly approach zero as T ^ 0 (pink curves). 


Types of a third group only appear at high temperatures 
(grey curves). Remarkably, these groups correspond ex¬ 
actly to the sets of primary A, secondary A, and non-A- 
types for each system, as enumerated using the analysis 
of instabilities, described above. The theory of A-types is 
thus consistent with observed results and suggests a new 
approach for analyzing thermal effects. 

One noteworthy feature of Fig. [^is the similarity be¬ 
tween FCC and HCP, in contrast to BCC. These gen¬ 
eral behaviors appear to depend primarily on the crys¬ 
tal structure rather than on bonding particulars. In¬ 
deed, preliminary investigations show that when atoms 
in BCC, FCC, and HCP crystals are perturbed from their 
equilibrium positions with a Gaussian noise whose mag¬ 
nitude scales with temperature (i.e., an Einstein model 
m), the distribution of topological types changes in a 
manner similar to that observed in Fig. 

A second notable feature is the total frequency of A- 
types in the liquid phases of the three systems. In liquid 
tungsten, the unique BCC-type accounts for less than 
0.05% of all atoms just above T^, where is the bulk 
melting temperature. In contrast, liquid copper consists 
of roughly 20% FCC-types, and liquid magnesium con¬ 
sists of roughly 30% HCP-types, just above Tm- These 
structural data are relevant in studying physical pro¬ 
cesses such as crystallization m and melting [22] . 

A third feature that stands out is the high fraction of 
A-types in the FCC and HCP samples at temperatures 
just below melting. In FCC, 94% of all atoms are classi¬ 
fied as having FCC local structure just below melting; in 
HCP this number is 96%. At 0.85Tm, these numbers are 
over 99% in both systems. This suggests a natural use 
of A-types for identifying structure in highly-perturbed 
atomistic systems, such as those at high temperatures. 

We next consider several applications of the proposed 
approach to illustrate some of its unique features. 


III. DEFECT VISUALIZATION 

As noted, the high frequencies of A-types in single crys¬ 
tals, even at extremely high temperatures, suggests their 
use for visualization of local structure in atomic systems. 
Figure shows thin cross-sections from an FCC alu¬ 
minum polycrystal prepared using MD [23]. The sam¬ 
ple was obtained by annealing a microstructure obtained 
through simulated grain-growth m, plastically deform¬ 
ing it at low temperature, and then thermalizing it at 
0.9Tm. In these figures, atoms that are FCC-types are not 
shown for clarity; among the remaining atoms, those that 
are HCP-types are shown in gold, and all other atoms are 
shown in dark blue. 

In Fig. I^a), grain boundaries can be identified as a 
network of non-FCC-type atoms. Vacancies (A) can be 
identified within the grain interiors as small clusters of 
non-FCC-types; as only a thin cross-section of the ma¬ 
terial is shown, not all atoms adjacent to these defects 
appear in the figure. A twin boundary (B) and stack- 
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FIG. 4. Frequencies of all A-types and all non-A-types that appear in the single-crystal (a) BCG tungsten, (b) FGG copper, and 
(c) HGP magnesium as a function of temperature upon heating from T = 0 to 150% of the bulk melting temperature. Blue 
curves indicate primary A-types, pink curves secondary A-types, and grey curves non-A-types. Thick curves indicate the sum 
of all frequencies of the corresponding color; note that there is only one primary A-type and no secondary A-types for BGG. 


ing fault (C) can be identified as one and two layers of 
gold-colored (HCP, non-FCC type) atoms, respectively. 
Figures l^b) and (c) show magnified images of a disloca¬ 
tion and a stacking fault inclined at a low angle relative 
to the cross-section plane. 

As noted earlier, individual Voronoi topologies can be¬ 
long to multiple families; we use the term indetermi¬ 
nate type to refer to such cases. This indeterminacy 
complicates the visualization procedure suggested here, 
as many types in Jmcp also belong to Tfcc • For this rea¬ 
son some atoms that belong to the twin boundaries and 
stacking faults are not seen in Figs. [^a) and (c). This 
shortcoming can be easily addressed within the topolog¬ 
ical framework, and is discussed in the supplementary 
material. 


The utility of the proposed topological framework for 
local structure classification and identification is useful 
for finite-temperature simulations of atomic systems con¬ 
taining defects. Of particular interest are the many phe¬ 
nomena which are only of interest at high temperature, 
such as dislocation climb |25], interface thermal fluctua¬ 
tion |26] , and defect kinetics under irradiation conditions 
EH- Conventional visualization methods require quench¬ 
ing or temporally averaging a sample prior to crystal de¬ 
fect analysis [2]. In general, we do not know whether 
this “tampering” with the data leads to significant dis¬ 
crepancies between the observed structures and the ac¬ 
tual finite-temperature ones. Remarkably, the proposed 
approach identifies all defects in Fig. [^and does not er¬ 
roneously identify any bulk atoms as being adjacent to 
defects, all at very high temperature and without quench¬ 
ing or time averaging. The topological approach thus 
provides a natural method for identifying and visualizing 
local structure that involves no ad hoc cut-off parameters 
and which is robust at high temperatures. This opens 
new opportunity for in situ structure analysis of atomic 
simulations at temperature, potentially identifying new 
high-temperature mechanisms and defect structures. 


IV. COMPARISON WITH OTHER METHODS 

Although a complete survey of existing methods for 
analyzing structure in high-temperature atomic systems 
is beyond the scope of this introductory paper, we briefly 
consider how visualization using A-types compares with 
three of the most frequently-used order parameters: cen- 
trosymmetry [28] , bond-angle analysis [29] , and adaptive 
common-neighbor analysis EIISq]. 

As a concrete example, we consider a stacking-fault 
tetrahedron m in a high-temperature FCC aluminum 
single crystal. A stacking-fault tetrahedron (SFT) is a 
three-dimensional defect consisting of four stacking faults 
that form the faces of a tetrahedron. The interior of an 
SFT is an FCC crystal; its edges are stair-rod disloca¬ 
tions [25]. Figure 1^ illustrates a cross-section through 
the center of an SFT and parallel to one of its four faces; 
the intersection of the SFT boundary with the view¬ 
ing plane is an equilateral triangle. This perfect SFT 
was constructed in FCC copper and then thermalized 
at 0.85Tm. Centrosymmetry, bond-angle, and adaptive 
common-neighbor analyses were all performed using the 
OVITO software package [32] , 

Figure]^ a) shows atoms colored using the centrosym¬ 
metry order parameter. In this coloring, atoms belong¬ 
ing to faces of the SFT have higher centrosymmetry val¬ 
ues than those in the FCC environment, as expected. 
Note, however, that many atoms inside and outside the 
SFT also have high centrosymmetry values. Such atomic 
environments do not, however, indicate crystal defects, 
but rather result from thermal fluctuations which locally 
distort the atomic environment. The inability of the 
centrosymmetry order parameter to distinguish struc¬ 
tural defects from thermal perturbations requires users 
to quench a system before analyzing it. 

Figures [^b) and (c) show atoms colored using bond- 
angle analysis and adaptive common-neighbor analysis, 
respectively, also at 0.85Tm. In these figures, many atoms 
belonging to the SFT faces are classified as having HCP 
local structure, as expected. However, both methods 
identify many atoms away from the stacking faults as 
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(b) Dislocation (c) Stacking fault 


FIG. 5. Polycrystalline aluminum at 938 K (0.9Tm); the width 
of each cross-section is 2 nm. Atoms that are FCC-types are 
not shown for clarity. Of the ones remaining, those that are 
HCP-types are shown in gold and all other atoms are shown 
in dark blue. Grain boundaries are seen as a network of non- 
FGG-types (dark blue and gold atoms). In cross-section (a), 
defects are labeled as follow: vacancies A, twin boundary B, 
and stacking fault C. Gross-sections (b) and (c) show magni- 
hed images of a dislocation and stacking fault. 


structural defects, despite the absence of other defects in 
the crystal. Moreover, application of bond-angle analysis 
incorrectly identifies many atoms in the bulk as having 
HCP local structure. Although the general shape of the 
SFT can be discerned in both figures, the details are am¬ 
biguous and automated location of the SFT in an atomic 
ensemble is difficult or impossible at the simulation tem¬ 
perature. 

These results are in contrast with the picture produced 
using Voronoi topology, illustrated in Figure |^d). The 
approach taken here provides the clearest representation 
of the SFT. In this case, every atom characterized as be¬ 
ing in an HCP environment is on an SFT face, without 
exception, despite the high temperature and the strain 
fields of the constituent partial dislocations. Moreover, 
all atoms not at the surface of the SFT are correctly 



(a) Gentrosymmetry 


(b) Bond-angle analysis 




(c) Adaptive GNA 


(d) Voronoi Topology 


FIG. 6. Gross-section of a stacking-fault tetrahedron in cop¬ 
per at 85% of its melting temperature, colored using several 
popular visualization approaches, and the proposed one. In 
(b) and (c), dark blue, yellow and red indicate atoms in FGG, 
HGP, and other local environments, respectively. In (d), dark 
blue, yellow, and red indicate atoms that are FGG-types, 
HGP- but not FGG-types, and all other types. 


identified as being FCC-type. Finally, atoms lying at the 
corners of the triangular cross-section through the SFT 
triangle are identified as neither stacking faults (HCP- 
type) nor bulk type, but as having a distinct local struc¬ 
ture; these are the dislocation cores. The sole weakness 
of this visualization procedure results from indeterminate 
types which belong to both J^fcc and Jfficp and whose 
neighborhoods are identified as FCC rather than HCP. 
This limitation can be addressed and is discussed in the 
supplementary material. 

This topological approach to structure visualization 
can also be applied to low-temperature systems such as 
those obtained through quenching (inherent structures); 
an example can be found in the supplementary material. 


V. GRAIN-BOUNDARY 
CHARACTERIZATION AND ANALYSIS 


The proposed framework also enables the analysis of 
structurally-complex systems in an automated manner. 
To illustrate this capability we analyze how a particular 
grain boundary transforms between a series of distinct 
structures as a result of absorbing point defects, as it 
may, for example, under irradiation conditions. 

In particular, we consider a 115 [001] (310) symmet¬ 
ric tilt boundary in a BCC tungsten bicrystal. This 
grain boundary exhibits three structurally distinct, sta¬ 
ble phases. We begin by characterizing these phases 
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FIG. 7. (a-c) Three structurally distinct stable phases of the 
E5 [001] (310) symmetric tilt boundary in BCC tungsten, ar¬ 
bitrarily labeled Phase I, Phase II, and Phase III. Atoms are 
colored according to topological type, each assigned a unique 
color; BCC-type atoms are not shown for clarity, (d) An 
image from a large simulation of the same boundary, initially 
constructed uniformly of the Phase I structure, and into which 
self-interstitial atoms have been randomly placed at a con¬ 
stant rate at 1500 K. Atoms are colored by topological type, 
as shown in (a-c); BCC-type atoms are not shown for clarity; 
all other atoms are shown in grey. 


using the Voronoi topologies of the constituent atoms; 
these phases are illustrated in Fig. [^a-c). Atoms are 
colored according to their topological type; BCC-type 
atoms are not shown for clarity. Phase I consists of three 
distinct topological types, colored in different shades of 
blue, Phase II consists of two distinct topological types, 
colored in shades of green, and Phase III consists of six 
topological types, colored in shades of red, orange, and 
yellow. 

We initialize the simulation by constructing a large 
bicrystal containing a E5 [001] (310) symmetric tilt 
boundary at 0 K and uniformly of the Phase I structure. 
The sample is then heated to 1500 K, and equilibrated at 
this temperature for 4 ns. Self-interstitial atoms are then 
inserted into random locations in the boundary plane at 
a rate of 1.62 atoms/nm^/ns. We analyze the result¬ 
ing grain boundary structure throughout the MD simu¬ 
lation using Voronoi topology. Figure [^d) shows a grain 
boundary with distinct domains of all three grain bound¬ 
ary phases, suggesting a phase transition driven by ab¬ 
sorption of self-interstitial atoms. 

To study the transformation of the grain boundary 


structure, we track the fraction of each phase present 
during the evolution. We begin by counting the number 
of atoms in the sample with topological types associated 
with each of the three phases. We next calculate the 
number of non-BCC-type atoms per nm^ in each of the 
three phases. Finally, we normalize the A-type counts for 
the three phases by dividing by the number of A-types 
per unit area and the total grain boundary area. 

Figure shows the fraction of the three phases over 
time, starting when the first self-interstitial atom is 
added to the grain boundary. During the initial 100 ps, 
there is a sharp decrease in Phase I, accompanied by a 
rapid growth of Phase III. After approximately 300 ps, 
the grain boundary structure settles into a pattern of in¬ 
creasing and decreasing Phase I, Phase II and Phase III 
fractions, all with the same period. The minimum in 
the Phase III fraction corresponds to the maximum in 
the Phase I fraction and the maximum in the Phase II 
fraction corresponds to minima in the Phase I and III 
fractions. The period is commensurate with the time 
required to add a full (310) plane of atoms to the sam¬ 
ple. At this temperature, the equilibrium grain boundary 
structure is dominated by Phase III, Phase I (the equi¬ 
librium structure at 0 K) never completely disappears, 
and Phase II is present only in a very small fraction of 
the grain boundary. 



FIG. 8. The fraction of the E5 [001] (310) symmetric tilt 
boundary in BCC tungsten occupied by its structurally dis¬ 
tinct, stable phases during the insertion of self-interstitial 
atoms at a hxed rate. 


This example illustrates the power of the topological 
approach for automating structure analysis in systems 
with complex defect structures and for determining de¬ 
fect statistics. We defer a more complete analysis of this 
example to a future study of grain boundary structure 
evolution during irradiation. 








































VI. DISORDERED STRUCTURES 

Finally, we consider how the proposed approach can 
be used to characterize disordered systems such as liquids 
and glasses. In contrast to conventional order parameters 
- which are typically useful for studying either ordered or 
disordered systems, but not both - the approach taken 
here can be applied effectively to all kinds of systems. 
As the topological type of each Voronoi cell provides a 
robust structural description of the local neighborhood 
of a particle, the distribution of topological types in a 
system can be interpreted as a statistical-topological de¬ 
scription of the system as a whole. This ability to char¬ 
acterize both ordered and disordered systems within the 
same framework is of particular importance for studying 
phase transitions between ordered and disordered phases, 
as well as between distinct disordered phases. 

Using MD, we simulate two disordered systems of cop¬ 
per atoms: a high-temperature liquid (HTL) equilibrated 
at roughly 1.85Tm, and a glass-forming liquid (GFL) 
obtained by undercooling the initial liquid to roughly 
0.75Tm; each system contained 1,372,000 atoms. The dis¬ 
tributions of topological types in the two systems enable 
us to describe structural features of the systems in a ro¬ 
bust and quantitative manner, and to observe structural 
differences between them. 



FIG. 9. Frequencies of the twenty-five most common topo¬ 
logical types in liquid copper at 1.85Tm, and corresponding 
frequencies in the glass-forming liquid copper at 0.75Tm; cir¬ 
cles indicate frequencies in quenched samples. 

Bars in Figure show frequencies of the fifty most 
common topological types in the HTL, and correspond¬ 
ing frequencies in the GFL. The most common types in 
the HTL do not occur nearly as frequently as the most 
common types in the GFL. In particular, the sum of fre¬ 
quencies of the five most common topological types in 
the GFL is greater than the sum of frequencies of the 
fifty most common types in the HTL. In this sense, the 
GFL is substantially more ordered than the HTL. Note 
that in an ideal gas, in which there is no interaction be¬ 


tween particles, the distribution of types is considerably 
less concentrated than in either of these systems [33] . 

Gircles in Figure indicate corresponding frequen¬ 
cies in systems obtained by quenching these two sys¬ 
tems. When quenched, frequencies of these common 
types change only modestly. This again indicates that 
the current method is relatively insensitive to thermal 
vibrations. 

Other means of quantifying disorder in atomic systems 
have been widely developed, and have been used to dis¬ 
tinguish distinct types of disordered systems. Structural 
correlation functions [T] and configurational entropy [34]- 
l36] are both non-local descriptors that have been used to 
study disordered systems. Other recent work has focused 
on local structural measures laiin]; our approach is in 
this vein. 


VII. CONCLUSIONS AND DISCUSSION 

We have introduced a mathematical approach towards 
classifying and identifying local structure of a particle 
within a system of particles. Applications highlight its 
utility in analyzing atomistic data sets, such as those pro¬ 
duced by molecular dynamic simulations. In particular, 
the theory of A-types enables identification and visual¬ 
ization of defects in ordered systems at high tempera¬ 
tures. This capability can play an important role for 
in situ study of high-temperature mechanisms currently 
inaccessible to current structure-identification methods. 
The proposed framework also enables a new approach 
for characterizing and identifying defects. This in turn 
allows for an automated approach for studying systems in 
which structural features evolve. Finally, Voronoi topol¬ 
ogy enables the characterization of disordered systems in 
a statistical manner, through the distribution of topolog¬ 
ical types. We have illustrated the potential of this ap¬ 
proach in distinguishing a high-temperature liquid from 
a glass-forming liquid. 

Any description of structure within a fixed distance of 
a particle will be unable to capture all long-range struc¬ 
tural features of a system. Figure provides clear ex¬ 
amples of this limitation, where an atom with HGP local 
environment might be part of a twin boundary, stack¬ 
ing fault, or other defect; further analysis is required to 
distinguish between these. The analysis of local struc¬ 
ture introduced here can be integrated into tools such as 
those developed in to automate long-range structural 
analysis. 

The authors have developed software to automate this 
analysis, and is available upon request. 


MATERIALS AND METHODS 

Deciding whether two Voronoi cells have the same 
topological type is equivalent to deciding whether two 
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planar graphs are identical, as the edge-boundary of ev¬ 
ery Voronoi cell is a planar graph. For each particle in 
a system we compute a “code” that records the graph 
structure of the edge-boundary network of its Voronoi 
cell. To do this, we first determine the Voronoi cell us¬ 
ing the Voro++ software package [38], which computes a 
list of faces, each represented as an ordered sequence of 
vertices. Next, we use a graph-tracing algorithm to com¬ 
pute a code for this planar graph. More specifically, the 
following algorithm of Weinberg m is followed: (a) An 
initial vertex is chosen and assigned the label 1. (b) An 
edge adjacent to that vertex is chosen and travel begins 
along that edge, (c) If an unlabeled vertex is reached, 
it is labeled with the next unused integer and we “turn 
right” and continue, (d) If a labeled vertex is reached 
after traveling along an untraversed edge, we return to 
the last vertex along the same edge but in the opposite 
direction, (e) If a labeled vertex is reached after traveling 
along an edge previously traversed in the opposite direc¬ 
tion, we “turn right” and continue; if that right-turn edge 
has also been traversed in that direction, we turn along 
the next right-turn edge available; if all outgoing edges 
have been traversed, we stop. At this point, each edge in 
the graph has been traversed once in each direction; the 
ordered list of the vertices visited is called a code. 

Codes are constructed for each choice of initial ver¬ 


tex and edge, and for each of two spatial orientations; 
all labels are cleared before producing each code. For a 
Voronoi cell with e edges, 4e codes are generated, each 
an ordered list of 2e integer labels. Each code completely 
describes the Voronoi cell topology, and so it is sufficient 
to only record one of them. A code for a typical Voronoi 
cell requires less than 100 bytes of storage. Additional 
details can be found in naiHisn]. 

Run-time. The use of Voronoi topology for structure 
identification is computationally efficient. In preliminary 
tests, the Voronoi topology of one millions atoms could 
be calculated on a desktop computer in less than one 
minute. By contrast, conventional approaches used in 
high-temperature structure analysis require that systems 
be quenched before visualization. A complete quench 
necessary to obtain the inherent structure can require 
several hours of computation for a system of comparable 
size. 
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SUPPLEMENTARY MATERIAL 


ENUMERATING PRIMARY 
AND SECONDARY TYPES 

In the paper we considered families of topological types 
IFx associated with particular structures. In this section 
we provide some additional detail regarding the determi¬ 
nation of these families, and report the numbers of types 
in several of them, as well as in the overlap between mul¬ 
tiple families. 

The Voronoi cell of BCC is topologically stable in the 
sense that infinitesimal perturbations of atomic coordi¬ 
nates will not change its topology. Therefore, J^bcc con¬ 
sists of a single primary type and no secondary types. 

In contrast, Voronoi cells of FCC and HCP are un¬ 
stable, and infinitesimal perturbations of atomic coor¬ 
dinates will change their topologies. The instability of 
these Voronoi cells can be detected in vertices that are 
incident with four edges; there are six such unstable ver¬ 
tices in FCC and HCP. As each unstable vertex can either 
remain unstable, resolve in one of 3 primary directions, 
or resolve in one of 4 secondary directions (see Fig. 3), 


we must check 8^ = 262,144 configurations that can re¬ 
sult from all infinitesimal perturbations. We compute 
the topology of each configuration using the algorithm 
described in the Materials and Methods section of the 
paper. Multiple configurations can result in the same 
Voronoi cell topology due to symmetries of the unper¬ 
turbed configuration. 

For FCC we find 44 primary types and 6250 secondary 
types; of the secondary types only 2771 have no unstable 
vertices. Figure illustrates several Voronoi cells ob¬ 
served in a finite-temperature FCC crystal; their topolo¬ 
gies are given by this enumeration technique. For HCP 
we find 66 primary types and 21,545 secondary types; of 
the secondary types only 9490 have no unstable vertices. 

While the determination of topological types associ¬ 
ated with a particular structure may require substantial 
computation, this needs only be done once per struc¬ 
ture. Lookup tables are then created and subsequently 
referenced when analyzing atomistic data sets. Data for 
several common crystal structures are available from the 
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FIG. SI. Examples of Voronoi cells in a perturbed FCC crys¬ 
tal. 


authors upon request. 

In the paper we noted that certain types can belong to 
several families (indeterminate types). More specifically, 
we note that the unique type belonging to Tbcc also be¬ 
longs to *Ffcc and Jmcp- Furthermore, J^fcc and Jmcp 
share 23 primary types in common and 1352 secondary 
types in common. 


RESOLVING INDETERMINATE TYPES 

In the paper we noted practical challenges created by 
the overlap of multiple families of types. In particular, 
when attempting to visualize defects in an FCC environ¬ 
ment, some atoms were mistakenly identified as belong¬ 
ing to the FCC bulk instead of to the HCP-like defect. 
Although a thorough analysis of this problem is beyond 
the scope of this paper, we briefly consider one approach 
towards resolving it. We leave a complete discussion of 
this topic for a future paper. 

Indeterminate types are Voronoi topologies that can 
be obtained through infinitesimal perturbations of mul¬ 
tiple distinct structures. Figure 2 in the paper shows a 
number of regions in configuration space incident with 
both xpcc and xhcp- One way of deciding whether a 
point in this region should be classified as FCC or HCP 
involves computing distances in this configuration space: 
points close to xpcc should be classified as FCC, while 
those close to xhcp should be classified as HCP. Points 
closer to xpcc than to either xpcc or xhcp should be 
classified as BCC. 

While theoretically appealing, there is no practical 
manner in which to calculate these distances because 
of the complicated topology of this configuration space. 
More specifically, the standard metric on cannot be 
used due to the actions of the rotation, renormalization, 
and permutation groups acting on it. 

A more practical approach involves perturbing con¬ 
figurations of particles as follows. If the Voronoi type 
of a particle is indeterminate, we randomly perturb the 
particle and its nearby neighbors; this corresponds to 
a small perturbation of x in configuration space. The 
Voronoi cell of the perturbed configuration is calculated; 
this procedure is repeated several times. If all pertur¬ 
bations result in indeterminate or HCP-types (the latter 
occurring at least once), then the particle is classified as 
being HCP-type. 

We consider a stacking-fault tetrahedron (SFT) in a 



FIG. S2. A high-temperature single crystal containing an 
SFT, visualized using (a) the method described in the paper 
and (b) the modified method proposed here. 


single copper crystal that was heated to 85% of its melt¬ 
ing temperature; this system was considered in “Compar¬ 
ison with Other Methods” section of the paper. Figure 


S2 shows the SFT colored by the approach considered in 


the paper and by the modification considered here. This 
modified approach provides a more robust visualization 
of the SFT than the approach suggested in the paper. 
We defer a complete discussion of indeterminate types 
and methods of resolving them to a future paper. 


QUANTITATIVE COMPARISON 
AT HIGH TEMPERATURES 


In the paper we have shown that Voronoi topol¬ 
ogy enables robust visualization of structure in high- 
temperature systems that cannot be obtained using con¬ 
ventional methods. Here we provide a direct and quan¬ 
titative comparison beyond visual inspection. We be¬ 
gin with a system containing 1,372,000 copper atoms 
organized in a perfect crystal, and heat the system to 
just below its melting point as described in the “Finite- 
Temperature Crystals” section of the paper. For each 
order-parameter considered in the paper, we calculate the 
frequency of atoms in this single crystal that are charac¬ 
terized as non-FCC-type as a function of temperature. 
Close inspection of the system shows that there are no 
defects. 

For each visualization method we considering the num¬ 
ber of atoms in this system characterized as non-FCC- 
type. We use OVITO to compute bond-angle analysis, 
adaptive common-neighbor analysis, and centrosymme- 
try values for each atom and at each temperature. Fig¬ 
ure shows the frequency of non-FCC-type atoms, as 
measured by each of the order-parameters, as a function 
of temperature. 

In contrast to the bond-angle and adaptive common- 
neighbor classifications provided by OVITO, centrosym- 
metry (CS) is reported as a scalar requiring choice of a 
cutoff; atoms with CS values above this cutoff are not 
considered FCC. To determine an appropriate cutoff we 
considered CS values of atoms adjacent to either a va¬ 
cancy or an interstitial in an otherwise perfect FCC crys- 
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FIG. S3. Frequency of atoms characterized as non-FCC-type 
in a copper crystal. 


tal. Atoms adjacent to a vacancy have a CS value of 6.20; 
interstitial atoms have a CS value of 13.88. In order that 
our choice of cutoff allows for the detection of vacancies, 
we choose a cutoff of 6. We note that cutoff values as low 
as 0.5 and 1 can be regularly found in the literature, and 
so our choice of cutoff is extremely conservative. 

At T = 1300 K, bond-angle analysis, adaptive 

common-neighbor analysis, and centrosymmetry incor¬ 
rectly identify over 36.9%, 55.7% and 13.5%, respectively, 
of the atoms as not belonging to the FCC bulk crystal. 
By contrast, Voronoi topology mistakenly identifies only 
0.53% as such atoms as not belonging to the FCC bulk 
crystal. 


QUENCHED SYSTEMS 

The topological approach is not limited to studying 
finite-temperature systems, and can also be applied to 
structures obtained through quenching high-temperature 
samples. Figure [S4| shows the same stacking-fault tetra¬ 
hedron (SFT) considered in the paper after quenching to 


0 K, colored using the same order-parameters considered 
in the paper. 

Although all methods successfully detect the pres¬ 
ence of the SFT in the quenched system, centrosymme¬ 
try, bond-angle analysis, and adaptive common-neighbor 
analysis also identify atoms outside the SFT as being 
non-FCC-type. These non-FCC-type identifications are 
the result of small local strains which are not structural 
defects. Visualization of the sample through Voronoi 
topology does not result in any such misidentideations, 
providing the most robust picture of structural defects. 



(a) Centrosymmetry 


(b) Bond-angle analysis 



(c) Adaptive CNA 


(d) Voronoi Topology 


FIG. S4. Cross-section of a stacking-fault tetrahedron in cop¬ 
per, after quenching from 85% of its melting temperature, 
colored using several popular visualization approaches, and 
Voronoi topology. In (b) and (c), dark blue, yellow and red 
indicate atoms in FCC, HCP, and other local environments, 
respectively. In (d), dark blue, yellow, and red indicate atoms 
that are FCC-types, HCP- but not FCC-types, and all other 
types, respectively. 






























